The Effect of Seasonal Variation and Hive Identity on Cold Tolerance in Tetragonula carbonaria

R code

The author acknowledges the assistance of AI tools in the preparation of this document (GitHub Inc. 2025; OpenAI 2025).

Loading Data

Code 
# Libraries
library(readxl)
library(dplyr)
library(knitr)
library(kableExtra)
library(moments) 
library(ggplot2)
suppressWarnings(library(plotly))
library(lme4)
library(performance)
library(sjPlot)
library(lmerTest)
library(DT)

# Excel data originating from the biodiversity_abundance data sheet
mydata <- read_excel("Data.xlsx")

# Convert categorical variables to factors
mydata$Month <- factor(mydata$Month)
mydata$Hive <- factor(mydata$Hive)
Code 
# Data table
datatable(
  mydata,
  options = list(scrollX = TRUE),
  rownames = FALSE
)
Table 1: Data for T. carbonaria exploring chill coma onset time (s) as a proxy for critical thermal minima (CTmin) across hive and month (n = 60 total observations).

Exploratory Data Analysis

Code 
# Summary statistics for chill coma onset time by hive and month
summary <- mydata %>%
  group_by(Hive, Month) %>%
  summarise(
    Mean = mean(`Coma time`),
    Median = median(`Coma time`),
    IQR = IQR(`Coma time`),
    SD = sd(`Coma time`),
    SE = SD / sqrt(n()),  # Standard error
    Skewness = skewness(`Coma time`),
    Kurtosis = kurtosis(`Coma time`),
    n = n(),
    .groups = "drop"
  )

# Summary statistics for chill coma onset time by month
overall_summary <- mydata %>%
  group_by(Month) %>%
  summarise(
    Hive = "Overall",
    Mean = mean(`Coma time`),
    Median = median(`Coma time`),
    IQR = IQR(`Coma time`),
    SD = sd(`Coma time`),
    SE = SD / sqrt(n()),  # Standard error
    Skewness = skewness(`Coma time`),
    Kurtosis = kurtosis(`Coma time`),
    n = n(),
    .groups = "drop"
  )

# Combine and arrange
full_summary <- bind_rows(summary, overall_summary) %>%
  arrange(Month, Hive)

# Render table
kable(full_summary, digits = 2)
Table 2: Summary statistics for chill coma onset time (s) for T. carbonaria used as a proxy for critical thermal minima (CTmin) across hive and month (n = 60 total observations).
Hive Month Mean Median IQR SD SE Skewness Kurtosis n
Hive 1 April 107.20 113.37 27.97 21.13 6.68 -0.07 2.07 10
Hive 2 April 122.22 119.36 14.73 27.66 8.75 1.25 4.75 10
Hive 3 April 108.52 111.88 20.63 19.00 6.01 -1.27 4.17 10
Overall April 112.65 114.30 25.55 23.15 4.23 0.69 5.61 30
Hive 1 May 141.68 145.20 30.71 28.17 8.91 -0.19 2.29 10
Hive 2 May 148.28 152.70 26.00 27.05 8.56 -1.50 4.77 10
Hive 3 May 189.27 164.82 85.41 74.69 23.62 0.59 2.21 10
Overall May 159.74 154.02 30.85 51.61 9.42 1.49 5.74 30
Code 
colour_scheme <- c("April" = "darkgreen", "May" = "lightgreen")
summary_df <- mydata %>%
  group_by(Month) %>%
  summarise(
    mean = mean(`Coma time`),
    se = sd(`Coma time`) / sqrt(n())
  )

# Column graph
graph <- ggplot(summary_df, aes(x = Month, y = mean, fill = Month)) +
  geom_col(width = 0.4, color = "black") +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2) +
  labs(
    y = "Chill Coma Onset Time (s)",
    x = "Month"
  ) +
  scale_fill_manual(values = colour_scheme) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05)), limits = c(0, NA)) +  # y-axis starts at zero
  theme_minimal(base_size = 14) +
  theme(
  panel.grid = element_blank(),
  axis.line.x = element_line(color = "black", linewidth = 0.7),
  axis.line.y = element_line(color = "black", linewidth = 0.7),
  axis.ticks = element_line(color = "black"),
  legend.position = "right"
)
ggplotly(graph)
Figure 1: Mean chill coma onset time (s) of T. carbonaria, used as a proxy for critical thermal minima (CTmin), compared between the warmer month of April and the cooler month of May. The error bars represent standard error. Each month includes 30 individuals (n=30).
Code 
# Boxplot
boxplot <- mydata %>%
    ggplot(aes(x = Hive, y = `Coma time`, fill = `Month`)) +
    geom_boxplot(position = position_dodge(width = 0.8)) +
    labs(
        x = "Hive ID",
        y = "Chill Coma Onset Time (s)",
    ) + scale_fill_manual(values = c("April" = "darkgreen", "May" = "lightgreen")) + theme(
  panel.grid = element_blank(),
  axis.line.x = element_line(color = "black", linewidth = 0.7),
  axis.line.y = element_line(color = "black", linewidth = 0.7),
  axis.ticks = element_line(color = "black"),
  legend.position = "right")
ggplotly(boxplot)  %>%
  layout(
    plot_bgcolor = "white",   # plot area background
    paper_bgcolor = "white"   # overall background
  )
Figure 2: Boxplot displaying the distribution of chill coma onset time (s) as a proxy for critical thermal minima (CTmin) in T. carbonaria across hives and months. The box plot is faceted by hive (Hive 1, 2 and 3) and coloured by month (April and May). Sample size is 10 individuals per hive (n = 60). The black dots represent outliers.
Code 
# Histogram
month_colors <- c("April" = "darkgreen", "May" = "lightgreen")
histogram <- mydata %>%
  ggplot(aes(x = `Coma time`, fill = Month)) +
  geom_histogram(color = "black", bins = 10, position = "identity", alpha = 0.7) +
  labs(x = "Chill Coma Onset Time (s)", y = "Frequency") +
  facet_wrap(~ Hive, scales = "free_x") + 
  scale_fill_manual(values = month_colors) +
  theme_minimal() +  theme(
  panel.grid = element_blank(),
  axis.line.x = element_line(color = "black", linewidth = 0.7),
  axis.line.y = element_line(color = "black", linewidth = 0.7),
  axis.ticks = element_line(color = "black"),
  legend.position = "right")

# Display the plot
histogram
Figure 3: Histogram displaying the distribution of chill coma onset times (s) as a proxy for critical thermal minima (CTmin) in T. carbonaria across hives and months. The histogram is faceted by Hive ID and coloured by month. Sample size is 10 individuals per hive (n = 60).

Assumptions

Code 
# Model assumptions
model <- lmer(`Coma time` ~ Month + (1 | Hive), data = mydata)
performance::check_model(model)

Figure 4: Graphical representations of the assumptions for the linear mixed model (n = 60) illustrating the distribution of residuals.
Code 
# Model assumptions for square-root transformation
mydata$Coma_time_sqrt <- sqrt(mydata$`Coma time`)
model_sqrt <- lmer(Coma_time_sqrt ~ Month + (1 | Hive), data = mydata)
performance::check_model(model_sqrt)

Figure 5: Graphical representations of the assumptions for the square-root transformed linear mixed model (n = 60) illustrating the distribution of residuals.

Model

Code 
# Model
model_lmerTest <- lmer(`Coma time` ~ Month + (1 | Hive), data = mydata)
tab_model(
  model_lmerTest,
  show.se = TRUE,
  show.stat = TRUE,
  show.ci = FALSE,
  show.re.var = TRUE,
  dv.labels = "Chill Coma Onset time (s)",
  string.pred = "Predictor",
  string.est = "Estimate"
)
Table 3: Summary of the linear mixed model predicting chill coma onset time (s) with random intercept for hive. Each month includes 30 individuals.
  Chill Coma Onset time (s)
Predictor Estimate std. Error Statistic p
(Intercept) 112.65 8.71 12.93 <0.001
Month [May] 47.09 10.16 4.63 <0.001
Random Effects
σ2 1549.42
τ00 Hive 72.73
ICC 0.04
N Hive 3
Observations 60
Marginal R2 / Conditional R2 0.258 / 0.291

References

GitHub Inc. 2025, GitHub Copilot [Code completion tool], GitHub Inc., https://github.com/copilot.

Native Bee Hives (n.d.) Tetragonula carbonaria [Photograph]. Available at: https://www.nativebeehives.com/wp-content/uploads/photo-gallery/nativebeedaisy1.jpg (Accessed: 27 May 2025).

OpenAI 2025, ChatGPT (Version 4), OpenAI, viewed 18 May 2025, https://openai.com/chatgpt.